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INTERPOLATION OF MSS DATA 

C. D. McGillem 

1. Introduction 

When continuous data Is discretized through sampling there is 
inevitably distortion that occurs when the data is reproduced. Further- 
more, there is no simple way to increase the scale of the sampled data 
set. A good example of a data set having this problem is the output of 
LiVe ERTS Mul t i spectra 1 Scanner System. This data set is sampled at a 
discrete set of spatial coordinates and is quantized in amplitude. When 
the data is reproduced for visual observation, there Is an intrinsic 
graininess in the picture due to the sampling. When an attempt is made 
to enlarge the picture, difficulty is encountered because data exists 
only at discrete spatial coordinates and there is no data between these 
points. One method of attacking this problem is to repeat each point a 
number of times thus generating a picture that is enlarged in proportion 
to the number of repetitions of the points. As can be imagined this 
type of enlargement does not lead to an improvement of the picture but 
rather serves to make some of the details more easily seen. Figure 1 
shows this type of enlargement of a portion of an ERTS frame taken 
from the Washington, D.C. region. Figure la is a reproduction of the 
data using each point only once. Figure lb is a A x A enlargement in 
which each point is repeated sixteen times. The graininess and artificial 
appearance of the enlarged picture are very evident on close inspection. 
From a distance the appearance is somewhat more acceptable but still 
has a very artificial look. 

An alternative method of producing an enlarged picture is to compute 
new values between tha original sample points by means of an appropriate 
kind of interpolation, it is the purpose of this paper to describe three 
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such interpolation procedures and to illustrate the result? obtained by 
applying them to the same area as shown in Figure 1. 

The mathematical bases of the interpolation procedures will be 
presented first and then the results of applying these procedures will 
be given. 

2, Mathematical Bases of Interpolation 

There are various rationales for selecting interpolating functions. 
However, none has been found that indicates the appropriate function to 
use for interpolating the ERTS MSS data so as to minimize any errors 
between the interpolated image and the true original image. There are 
reasons to believe that such a rationale may exist and when it has been 
properly defined a new kind of interpolation may be appropriate. In 
the meantime three more or less conventional procedures will be used: 
polynomial interpolation, trigonometric interpolation and sine function 
interpolation. The first two methods involve passing an appropriate 
curve through the data points surrounding the range of coordinates where 
the interpolation Is to be performed and then computing the interpolated 
values from the curve. The third method consists of taking a weighted 
sum of all data points in the set to compute the interpolated values. 

The weighting function is the sine function and this type of interpolation 
is exact for samples taken froi a continuous function having no spatial 
frequency components higher than one half the sampling frequency. These 
procedures will be considered individually. 
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(a) Original Data Points 



(b) Repeated Points 4x4 
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Washington D.C. area Run 72041900 10/11/72 Ch 3, Lin 
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3. Polynomial Interpolation (PQLYINT) 

By passing a nth order polynomial through n + 1 points it is 
possible to compute interpolated values from the resulting 
polynomial. There are several equivalent ways of looking at this 
type of data modification. Among the most useful are: computing 

interpolated values as a weighted sum of surrounding points; 
two dimensional convolution; and processing in the frequency domain. 
These various approaches will be illustrated In the following 
example of interpolating three intermediate points between equally 
spaced samples. 

The technique of achieving two dimensional interpolation 
will be to carryout a sequential operation: interpolate in 

x-direction first and then interpolate the modified data in the 
y-direction. This process is equivalent to assuming that the 
interpolation surface can be represented as the product of functions 
involving only one coordinate, ?,e., if the two-dimensional 
interpolating surface is g(x,y) then it is assumed 

g(x»y) = g^x) g 2 (y) (i) 

The proper weighting for polynomial interpolation can be 
computed from the Lagrarge interpolation formula as follows: 


f(x) = L f (X ) + L, f(X.) + ... L f(Xj 
ooil n n 

where 

f(X^) = kth sample value 


n-1 


/ n-1 

= II 

(X-Xj) / 

n 

j=o 

j=0 

jVk 

/ 

j^k 


( 2 ) 


( 3 ) 
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For a cubic polynomial four sample values are required and 

the coefficients can be written as follows. (Note the change in 

indices so that x varies from x to x, in the central portion of 

o I 

the interpolation region.) 

(x-x ) (x-x. ) (x-x.) 

I „ o , . * , 

-1 “ (x_,- x o ) (x_ 1 “Xj ) (x_ r x 2 ) 

(x-x_ 1 ) (x-x 1 ) (x-x 2 ) 

L o ' lv x -i H V x i )( W 


(x-x , ) (x-x ) (x-x„) 

1 _ ^ ^ r . 

L 1 “ (x ] -x_ ] ) (x 1 -x o ) (x 1 -x 2 ) 

(x-x^) (x-X Q ) (x-x^) 

L 2 ” (x 2 -x_,) (x 2 -x q )'(x7 
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Considerable simplification results when the samples are 
equally spaced and when the coordinate x is expressed as a fraction 
of the sampling interval. For this case letting u = x/Ax the 
interpolation equation becomes. 

f(u) - - |u(l-u) (2-u)f(-l) + I(l+u){l-,:.'2-u)f(o) 

+ l(l+u)(u)(2-u)f(l) - £-(l+u)(i>; - ••u) f (2) (5) 

For a specific interpolation interval xhe coefficients can 

be evaluated and a specific equation determined. As an example 

consider the Interpolation of three equally spaced values between 

1 1 3 

the original values. For this case u = p j, ^ and the coefficients 


are 
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u 

L -1 

L o 

L 1 

L 2 

1A 

-7/128 

105/128 

35/128 

- 5/128 

1/2 

-8/128 

72/128 

72/128 

-8/128 

3A 

- 5/128 

35/128 

105/128 

-7/128 


The resulting equations are 
f(0) = f (0) 

f<T> “ 


-7'128f(-l) + -j|f f(o) + i^- f(l) 


f(i-) 


f (f) 


-8/l23f(-i) + f (o) + f ( 1 ) 


TTS" 

35 

T2F 


-5/128f(-l) + f (0) + f (1) 


riff f(2) 

TW f(2) 
tts f(2) 


( 6 ) 


This set of equations can be written in matrix form as 


f(0) 



= 





_ 

_ 


0 

1 

0 

0 

-7 

105 

35 

-5 

TTS 

ITS 

ITS 

TTS 

-8 

72 

72 

-8 

TTS 

nS 

T 2 F 

T 2 K 

-5 

TTS 

vh 

105 

TTS 

-7 

TY8" 



f(-l) 


f (0) 


f(D 


f( 2 ) 

_ 

— r 


(7) 


f = A f 
—11 n 


The process of interpolating can also be interpreted as discrete 
convolution. By augmenting the original time series with zeros 
at points where interpolated values will occur the interpolated 
time series can be written as 


«u aug) * thl 


( 8 ) 


(• • * f_2 0»0» f_ i ^,0,0,0,fQ,0,0,0,fj,0,0,0,f2>...)'" 

(0,-5,-8,-7,0,35,72,105,128,105,72,35,0,-7,-8,-5) x ^ 


( 9 ) 


This representation will be considered in more detail later when 
the three interpolation schemes are compared. 
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An equivalent process can be carried out in the frequency domain. 
Taking the discrete Fourier transform of the above equation we have 

i, N -i : )fr nk 

f< Tfr> ■ vw - 5 00, 

but since f(nAT/A) - 0 for all n not a multiple of A this can be 
simpl i f ied to 

N -1 ~J 2lf nk 

F D {f uaug } = E f(nT) e N k^O, 1 . . L/N- 1 (11) 

Since F(k^) is periodic in k with period of N it is evident that 
the above transform is just the original spectrum with three 
replications following it. 

In order to carry out the convolution by multiplication of 

the DFT it is necessary that the functions have the same number 

of points and that a sufficient number of zeros be added in the 

time domain to prevent alising. For the present instance this 

means adding 16 zeros after the data points. To keep the total 

number of points a power of two the number of data points should 
k 

be set equal to 2 - 16. For example if k = 7 the number of data 
points is 128 - 16 ■ 112 and the convolution is 
frp* f i * * ■ f 1 27 » 0 « « • 0) * (h Q ,h. ..hj,_,0,0. .0) 

The processing operation would then be 

{f inter } = {F(k) H(k)} (l2) 

where 3^ 1 is the inverse discrete Fourier transform operator and 
F(k) and H{k)are the discrete Fourier transform of f f { n ) } and 
(h(n)} respectively. These operations can be carried out using 
the fast Fourier transform algorithm. 

If the interpolation is to be done one line at a time it 
appears that matrix multiplication would be more rapid than using 
the FFT. However, if a number of lines are to be processed 
simultaneously the transform method may reduce the computation time. 


Trigonometric Polynomial Interpolation (TRIG I NT) 

By passing a nth order trigonometric polynomial through 
n + 1 data points an interpolation can be performed by computing 
intermediate values from the polynomial. This process is 
equivalent to computing the Fourier series expansion of a function 
having sample values corresponding to the data subset. One reason 
for using trigonometric interpolation is that the errors are 
uniformly distributed over the interpolation interval as compared 
to polynomial interpolation in which the errors are much greater 
near the end points.^ In the present case this is not a major 
consideration since the Interpolated values all lie in the center 
interval of a multi-interval data set. However, because it 
does utilize different basis functions for performing the inter- 
polation it was selected as one of the methods to be studied. 

The rate at which the coefficients of a Fourier series 
approximation of a particular function decrease is determined by 
the smoothness of the function. If the function is continuous 
and has a continuous first derivative, the coefficients decrease 

at least as rapidly as -L- where k is the order of the harmonic. 

k 

One way of obtaining these conditions for an arbitrary function 
represented by N sample values is as follows. 


(i ) Choose the interval (0,N) as one half the period 
of a period function. 


(ii) 


Subtract away from the data the linear trend from 
the first to the last point, i.e., form a new data 
set 


= y jL 


TT 


k + y 


( 13 ) 


This will made z = z„ - 0 and thereby make the periodic function 
o N 

cont inuous. 

(iii) Reflect the data set {z^} around the origin as an odd 


s 


function. This will cause the periodic function to have a continuous 
derivative at z u and z 


•N 


-N* 


(iv) Expand this new data set in a sine series according 
to the formulas. 


z(t) 

N- 1 

E 

k-1 

, . irk t 

b. sin r- 
k N 

04) 

b. 

2 

N ;' , ,r k 

y z 5 n ____ 

05) 

k 

N 

. n N 

n=d 

Both b and 
o 

b„ are 
N 

indeterminant since 

the sine functions are 

for both of 

these 

cases. 



Finally the interpolated values z(t) are computed and the 
linear trend of the data added in again giving 

y (t) = z(t) + — - t + y Q 

If the interpolation is restricted to the center Interval 

N 

then we can replace t by u ® t - — and allow u to vary from 0 
to 1. Thus, the equation for interpolation becomes 


(16) 


y (u) - z(y + u) + (' N 


y.. - y N 

N Y °) ( u + -£) + y Q 


07) 


As an example consider a six point interpolation for which N = 5 

Jj 

b, = -jr E z si n-^p- 
k 5 n-l n 5 

b 1 = O.A[.58779z ] + .95106z 2 - .58779z 3 - .95!06 Zi< ] 

b 2 » 0.A[.95106 Z] + .58779 z 2 - .587792, - .951 

b 3 - 0.4[.95106z 1 - .58779h 2 - .58779z 3 + .95106z A l 

b 4 = 587792, - .95106 z 2 + .95106z 3 - .58779z ij ] 

z(u) = b, sin ~ (u+2.5) + b 2 sin (u+2.5) + b^ sin •— (u+2.5) 

+ b^ sin (u+2.5) 

ir , . 2 tt , 3mr , . 4ir \ 

= b, cos — u - b sin u - b_ cos u + b. sin — u 
1 5 c 5 3 5 4 5 


08) 


09) 


(20) 


10 


N-l 

z(u) a E b sin 
k«l k 


irk (u+N/2) 
~“N 


( 21 ) 


N "' V 2 _ .."kn flk (u+N/2) 

t E . N 2 n 5ln T 5,n~ 

k=l n a l 


( 22 ) 


This can be put Into the weighting function form by combining all 
the coefficients that multiply the same sample value. The weighting 
function form Is as follows. 


P(u) s n 


N-l 

E 

k=1 


z n s?n IT sin 


irk (u+N/2) 
N 


n < u < n + 1 


(23) 


5. Sine Function Interpolation (SINCINT) 

When a signal is bandlimited to frequencies no higher than 
one naif of the sampling frequency it is possible to exactly 
reconstruct the signal from its samples. The expression Is 

tw 

x(t) = L’ x(k) sine { t- k) (2*0 

k=- £ " 

where sine t is defined as slnntArt. The above operation is 
clearly the convolution of sine t with samples of x(t) as discussed 
earl ter. 

When there are actually higher frequency components present 
this reconstruction leads to the bandlimited function that most 
closely fits the original data set. 

One problem that arises Immediately in using this Interpolation 
procedure Is the requirement for a very large data set. The 
reason for requiring a large number of points is the very slow 
rate at which the interpolating coefficients fall off away from 
the point being interpolated. The large width of the Interpolating 
pulse means that it is possible to get undesirable edge effects 
that will persist through the interpolated picture if there is a 
sharp discontinuity at the start of the image. These difficulties 
can be essentially eliminated by using either a Fourier series 
expansion of the data set or the discrete Fourier transform which 
is virtually the same thing. Using the latter approach the 
interpolation is accomplished as follows: 

(i) The DFT of the data set is calculated. The 
highest frequency component will be one half 
the sampling frequency. 
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(i!) Under the assumption that the data is bandl I mi ted 
the spectrum is extended to any arbitrary higher 
frequency by adding zeros for the higher frequency 
samples. 

(Hi) The augmented (with zeros) spectrum ts now inverted 
to give the original function but with data points 
spaced by an amount determined by the number of 
zeros added to the original spectrum. If there W'Xo 
originally N samples and kN zeros are added to the 
spectrum then there will be (k+l)N points in the 
reconstructed signal. Their the scaling will be 
k+1. 

This operation Is carried out most rapidly by means of the 
fast Fourier transform. In order to use the FFT the number of 
original data points must be an integral power of two and the 
scale factor k+1 must also be a power of two. 

If there is concern about the presence of large discontinuities 
in the data that may lead to anticipatory or trailing oscillations, 
it is possible to modify the Interpolation process by multiplying 
the spectrum by an appropriate Hamming window function that will 
taper the high-frequency response so as to reduce such oscillations. 
This capability is included in the interpolation routine and is 
provided as an option to the user. 
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6. Interpolation as a Convolution Operation 

Interpolation can be viewed as the result of convolving samples 
of a signal with an appropriate interpolating pulse. Comparison 
of the interpolating pulse shapes corresponding to different 
interpolation schemes provides considerable insight into their 
relative performances. The basis of the interpolation-convolution 
analogy is the fact that both may be considered as weighted sums 
of the original data points. Tor interpolation we have 

00 

f(u) = £ f(k) c. (u) ( 25 ) 

k““«° K 


= ...c_ 1 (u)f (-1 )+c q (u) f (o)+c ] (u)f (1)+... ( 26 ) 

where c^(u) is the interpolation coefficient corresponding to the 
point f(k). For example, in the case cubic polynomial interpolation 
there are only four non-zero coefficients, namely 


C_ 1 (u) 

= 

- g- u(l-u) (2-u) 


C o (u) 


i(l+u)(1-u)(2-u) 


C, (u) 

= 

j u(1+u) (2-u) 


C 2 (u) 

= 

- g- u(l+u) ( 1 — u) 

(27) 

convolution operation 

can be expressed as 


00 

f(u) - { £ 

f(u) 

6 (u-k) ) * p (u) } 

(28) 


K=-« 

“ Z f(k) p(u-k) (29) 

= . ..p(u+i ) f (-1 ) + p(u)f(o) + p(u-l)f (1)+... ( 30 ) 




J 


\b 


where p(u) is the interpolating pulse. From (26) and ( 30) we obtain 
the following correspondence between the Interpolating pulse and the 
weighting function. 


p(u-k) 

= t,. (u) 

r. 

0 < u < l 

p(u) 

«= c^u+k) 

0 < u+k < 1 
-k < u < 1-k 


For the cubic interpolation coefficients we obtain 
p(u) - c^(u+2) = - ^(u+2) (u+3) (~ l~u) 

- c ^ (u+1 ) = y(u+l ) (u+2) (1~u) 

= c (u) = --(l+u) (1-u) (2-u) 

O Z 

= c_ 1 (u-1) = - ^(u-l)(2-u)(3-u) 


-2 < u < -1 

-1 < u < 0 

0 < u < 1 

1 < u < 2 


Substituting u = -u into the above equations shows that p(u) is an 
even function of u and so only the positive values need be calculated. 

Similar calculations can be carried out for the other interpolation 
methods and the convolution operator determined. Figure 2 shows the 
convolution pulses for the three kinds of interpolation considered here: 
polynomial, trigonometric and sine function. 

There is clearly much similarity among the three Interpolation 
operations. The most noticeable difference is the extended nature of 
the sine function compared to the other two. Because of this it is 
possible to get significant edge effects if there is a large dis- 
continuity in amplitude at the edge of the image, it is evident that 
all of the interpolation schemes will give some overshoot when a step 
discontinuity in amplitude is encountered. However this has not been 
found noticeable in the processed images — possibly because of the 
limited number of gray shades that can be viewed on the display. 
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7. Examples of Interpolated Imagery 

A number of examples of interpolated F.RTS imagery are shown in 
the accompanying figures. Using POLYINT, the area in the vicinity of 
the Pentagon Building is shown with magn i f i cat ion of 4 x 4, 8x3 and 
32 x 32 in Figure 3. The smooth transition between points of greatly 
differing contrast is evident. No overshoot or ghosting is evident 
in these images. 

Figure A shows the Pentagon interpolated with SINCINT. In Figure 
4(a) the image Is displayed with 16 gray levels and is seen to be very 
similar to the comparable image obtained using POLYINT and shown in 
Figure 3(c). Figure Mb) and Me) show the effect of reducing the 
number of gray levels displayed and illustrate the type of contouring 
that can be obtained in this manner. 

Figure 5 shows the same area interpolated using TRIGINT. However, 
in this case a different magnification is used in the vertical and 
horizontal directions to correct for the scale differences of ERTS 
images in these directions. The scaling ratio used is 38:27 which 
is quite close to the correct value of 79:56 specified in the ERTS Data 
Users Handbook, Figure 5(a) shows 16 gray levels which Figure 5(b) 
and 5(c) are binary images with the thresholds set at the Mh and 3rd 
gray levels respectively of the original image. 

Figure 6 shows a comparison of TRIGINT and SINCINT using a different 
ERTS frame. Figure 6(a) is the original data set while 6(b), 6(c) and 
6(d) are A x 4 magnifications using repeated points, TRIGINT and SINCINT 
respectively. Figure 4(d) has a somewhat different appearance than 4(c). 
The texture appears to be mottled somewhat. This is a general charac- 
teristic of images interpolated by SINCINT and has not been satisfactorily 
explained as yet. One possibility is that the SINCINT process is pro- 


viding a type of enhancement in which adjacent points are more completed 
separated than in the other kinds of interpolation. 

It can be concluded from examination of the above examples of 
interpolated imagery that the scene appears to be enlarged without 
introducing any major changes in its appearance. In many cases details 
of the scene are more easily discerned. 
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(a) POLYINT W x Lin 1129-1257 Col 1217-1 3^5 



(b) POLYINT 8V x 8 h Lin 11^8-1212 Col 1213-1277 



(c) POLYINT 32V x 32H Lin 1172-1188 Col 12^2-1258 
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Fig. 3 Pentagon 

Run 720^1900 10/11/72 Ch 3 


(a) 


SI NCI NT 32 V x 32H 



(b) SINCINT 32 V x 32H k Gray Levels 



(c) SINCINT 32V x 32H 2 Gray Levels 

«>! : i.VAI. I’AOK It 
'it 1 1 1. H’ ' ‘I ■ 1TV 


Fig. I 4 Pentagon 

Run 720141900 10/11/72 Ch 3 
Lin 1172-1182 Col 121*2-1258 
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(a) TRIGINT 32 V x 27H 16 Gray Levels 



(b) TRIGINT 32V x 2 /H Threshold 4th Level 



(c) TRIGINT 32 V x 27H Threshold 3rd Level 


oruus.m. . -• > ; 
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Fig. 5 Pentagon 

Run 72041900 10/11/72 Ch 3 
Lin 1172-85 Col 1242-1260 



(a) Original Data 


(b) Repeated Points it x it 


(c) TRIGINT x it 


SI NCI NT It x ^ 


Fig. 6 Calument Area 
Line 1 326- 1^53 


Run 72059900 Chan. 2 
Col 1 360- 1 ii87 
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